%{
                                                          vimArea.m V1
                                                     Written by Rajeev Kant
                                                        Sept. 4th, 2020
                                    Coulombe Lab for Cardiovascular Regenerative Engineering
                            Center for Biomedical Engineering, Brown Univeristy, Providence, RI, 02903

This function takes in vimentin stains and creates masks based on size to
differentiate between single cells and networks, then calculates area.

The user may need to edit the following lines: 69 (particle removal
algorithm), 90-91 (isolated and network mask thresholds)

Inputs:
1. region: IF or R image of vimentin stain
2. nuclei: binary nuclei image from analyzeNuclei.m to fill in holes in
Vimentin mask, helps mitigate network breakage
3. ROIpos: Position information of subROI being run determined in getROIs.m
4. path: file path
5. regName: Character array for region ('IF' or 'RM') for naming

Outputs:
1. summary: Summary values for output
2. Saved images of single/network masks


Original publication for citation: 
Kant, R. J., Bare, C. F., and Coulombe, K. L. K. "Tissues with 
patterned vessels or protein release induce vascular chemotaxis
in an in vitro platform", Tissue Eng. Part A, 2021.
doi:10.1089/ten.TEA.2020.0269

%}
function [summary] = vimArea(region, nuclei, ROIpos, path, regName)
%% Define subROIs for Vimentin stains

% Set predefined ROI positions based on input image and ROIpos
subR1 = imcrop(region,'Position', ROIpos(1,:));
subR2 = imcrop(region,'Position', ROIpos(2,:));
subR3 = imcrop(region,'Position', ROIpos(3,:));
subR4 = imcrop(region,'Position', ROIpos(4,:));
subR5 = imcrop(region,'Position', ROIpos(5,:));
subR6 = imcrop(region,'Position', ROIpos(6,:));
vimROIs = cat(3, subR1, subR2, subR3, subR4, subR5, subR6);
subROI_px = ROIpos(1,3)^2; % subROI area size - subR# size is actually 1 row and column greater than specified ROIpos

%% Binarize and remove small particles based on size exclusion
% Alter bwareaopen and numel() multiplier as needed
% General procedure is the same as analyzeNuclei.m but optimized
% for Vimentin stains and without segmenation

% Initialize
masks = zeros(size(subR1)); % Total ROI binary masks
maskSing = zeros(size(subR1)); % Single cell masks
maskNet = zeros(size(subR1)); % Network masks
AROI = zeros (1,6);
Anet = zeros (1,6);
Asing = zeros (1,6);

for i = 1:6
    % Create two separate binary masks at different sensitivities to include full extent of network
    % and minimize breakage
    bin1 = imbinarize(vimROIs(:,:,i));
    s = graythresh(vimROIs(:,:,i));
    bin2 = imbinarize(vimROIs(:,:,i),'adaptive', 'Sensitivity', s);
    binary_mask = bin1 + bin2; % Combine binary results
    
    for j = 1:2
        % Particle removal process
        binary_label = bwlabel(binary_mask,4);% Label each identified structure in the binary mask
        for labels = 1:max(binary_label,[],'all') % Remove all particles less than .00005 * size
            count = numel(find(binary_label == labels));
            if (count < 0.00005*numel(vimROIs(:,:,i)))
                binary_label(binary_label == labels) = 0;
            end
        end
        binary_mask = (binary_label ~=0);
        binary_mask = bwareaopen (binary_mask, 50); % Secondary particle removal
        
        if j < 2 % Optional - after first particle removal, run binarization on "leftover" mask to repair overthresholding breakages
            X = vimROIs(:,:,i) - uint8(255*(binary_mask)); % Calculate difference between original grayscale and final binarization
            %figure; imshow(X);
            binX = imbinarize(X); % Rebinarize
            %figure; imshowpair(binary_mask, binX, 'montage'); % Compare binary mask
           binary_mask = bwmorph(binary_mask + binX, 'thin', 1);
        end
    end
    
    binary_mask = binary_mask + nuclei(:,:,i); % Combine filtered Vimentin stain with binarized nuclei
    binary_mask = imbinarize(binary_mask);  % Reset values to true binary
    
    
    % Create single cell and network masks
    binSing = xor(bwareaopen(binary_mask,0),  bwareaopen(binary_mask, 500));
    binNet = bwareaopen(binary_mask, 500);
    
    AROI(1,i) = sum(binary_mask, 'all'); % Calculate area of binarized image
    Asing(1,i) = sum(binSing, 'all'); % Calculate area of binarized image
    Anet(1,i) = sum(binNet, 'all'); % Calculate area of binarized image
    
    % Generate and save false-colored image of masks
    figure; imshowpair(binSing, binNet, 'falsecolor');
    iptsetpref('ImshowBorder','tight');
    saveas(gcf,fullfile(path, strcat(regName,'_Vim',num2str(i),'.tif')),'tiff');
    
    % Add binary images and masks to their own arrays
    if i == 1
        masks = cat(3, binary_mask);
        maskSing = cat(3, binSing);
        maskNet = cat(3, binNet);
    else
        masks = cat(3,masks, binary_mask);
        maskSing = cat(3, maskSing, binSing);
        maskNet = cat(3, maskNet, binNet);
    end
    
end
close all;

%% Calculate summary statistics for output
% Normalized by analyzed area: Area (px^2)/(px^2)

% Percent area  relative to ROI size
A1 = AROI(1) / subROI_px;
A2 = AROI(2) / subROI_px;
A3 = AROI(3) / subROI_px;
A4 = AROI(4) / subROI_px;
A5 = AROI(5) / subROI_px;
A6 = AROI(6) / subROI_px;

% Relative area of cells that is from isolated cells
As1 = Asing(1) / AROI(1);
As2 = Asing(2) / AROI(2);
As3 = Asing(3) / AROI(3);
As4 = Asing(4) / AROI(4);
As5 = Asing(5) / AROI(5);
As6 = Asing(6) / AROI(6);

% Relative area of cells that is from networked cells
An1 = Anet(1) / AROI(1);
An2 = Anet(2) / AROI(2);
An3 = Anet(3) / AROI(3);
An4 = Anet(4) / AROI(4);
An5 = Anet(5) / AROI(5);
An6 = Anet(6) / AROI(6);

% Total, Inner, and Outer averages
scaleArea_Inner = mean(AROI(1:3)) / subROI_px;
scaleArea_Outer = mean(AROI(4:6)) / subROI_px;
scaleArea = mean(AROI(1,1:6)) / subROI_px;

scaleSingle_Inner = mean(Asing(1:3)) / mean(AROI(1:3));
scaleSingle_Outer = mean(Asing(4:6)) / mean(AROI(4:6));
scaleSingle = mean(Asing(1,1:6)) / mean(AROI(1:6));

scaleNet_Inner = mean(Anet(1:3)) / mean(AROI(1:3));
scaleNet_Outer = mean(Anet(4:6)) / mean(AROI(4:6));
scaleNet = mean(Anet(1,1:6)) / mean(AROI(1:6));

NSRatio = scaleNet / scaleSingle;
NS_Inner = scaleNet_Inner / scaleSingle_Inner;
NS_Outer = scaleNet_Outer / scaleSingle_Outer;

summary = [scaleArea, scaleArea_Inner, scaleArea_Outer, A1, A2, A3, A4, A5, A6, NSRatio;...
    scaleSingle, scaleSingle_Inner, scaleSingle_Outer, As1, As2, As3, As4, As5, As6, NS_Inner;...
    scaleNet, scaleNet_Inner, scaleNet_Outer, An1, An2, An3, An4, An5, An6, NS_Outer];
end % End of vimArea

